/*******************************************************************************

                                APODFO

  Title:        Apodize a file of fobs single-particle data and 
		prepare Wilson-style plots.

  Author:       Hanna Szoke

  Date:         6/4/01 - adapted from EDEN for SPEDEN. 

*******************************************************************************/

#include "util.h"
#include "cellparams.h"
#include "dims.h"	
#include "apodize.h"
#include "sp.h"

static	real	*Fap ;
static	real	*sigap ;
static  real	*Int ;			/* full set of intensities, (amp^2) */
static  real	*ksq ;			/* full set of (1/d^2)  */
static  real	*sigsqI ;
static  real	*sigsq ;
static  real	*sratio ;
static 	real	asig, bsig ;		/* coefs for apodization of sigmas */

static	void	analyze_sigsq() ;
static	int	apodfo_startup() ;
static	void	apply_apod_fo() ;	
static	void	prepare_fo() ;		/* make weighted F^2, sig^2 */

void	apodfo_main(int argc, char *argv[])
{
     void	writeFolist(char *) ;
     void	write_HDFfobs(char *) ;

	/* in apodutil.c:	*/

     void	allocate_Bin_arrays(int) ;
     void	do_binning(int, real *, real *, real *) ;	
     real	get_delta(real) ;
     void	init_apod_arrays(int, real *, real *, real *) ;
     real	linfit(int, real *, real *) ;		
     void	make_wilson(real *) ;		
     void	make_selected_wilson(char *) ;	
     void       scalefobs(real *, real) ;
     real	scale_fo(char *, char *) ;
     void	show_wilson_plots() ;

     char	out_filename[MAXSTRING] ; /* name of output file for x/y data */
     char	name_fo[MAXSTRING], name_fc[MAXSTRING] ;
     real	delta1 ;
     real	dstdev1 ;
     real	B1 ;
     real	factor ;

  /*******************************************************
  Check for right ballpark; identify fobs file name.
  *******************************************************/

     if (argc > optind+3) 	/* too many arguments */
	  ballpark(caller) ; 
     else if (argc == optind+3)
          sprintf(sf_filename, "%s", argv[optind+2]) ;
     else {
	  prompt("\nPlease enter the name of an fobs file: ") ;
          sscanf(terminp, "%s", sf_filename) ;
	  sprintf(message, "Running Apodfo on fobs file %s", sf_filename) ;
	  printTwice(message) ;
     }

     strip_path(sf_filename, pwd_name) ;

     hello(caller) ;

  /*******************************************************
  Read the input fobs file using standard read function.
  Set up arrays in (h,k,l) reciprocal lattice triplets,
  Prepare a list of non-zero intensities, sigma squared*
  intensities, and nusq's.  
  *******************************************************/

     Nbin = apodfo_startup() ; 
     
     init_apod_arrays(Nfob, Int, ksq, sigsqI) ;
     prepare_fo() ;

  /*******************************************************
  Bin the information, gathering logs of average binned 
  intensities & weights.  Determine coefficients (a, b) 
  for fitting sigsq as a linear function of Int.
  *******************************************************/

     allocate_Bin_arrays(Nbin) ;
     do_binning(Nfob, Int, ksq, sigsqI) ;

     if (useSig)
         analyze_sigsq() ;

  /*******************************************************
  Do a linear fit to establish slope and intercept.  
  Prepare for display the original binned log(I^2) and 
  the linearly fitted version.
  *******************************************************/

     B1 = linfit(Nbin, w_lnIsq, apodwt) ; 
     delta1 = get_delta(B1) ; 
     make_wilson(&dstdev1) ;

     show_wilson_plots() ;
     delta = delta1 ;

  /*******************************************************
  Replace Fap and sigap by apodized version, insofar as 
  delta > 0.
  *******************************************************/

     if (delta > 0) {
          apply_apod_fo() ;

          init_apod_arrays(Nfob, Int, ksq, sigsqI) ;
          prepare_fo() ;

          do_binning(Nfob, Int, ksq, sigsqI) ;
   	  sprintf(name_fo, "%s_wil", pwd_name) ;
	  make_selected_wilson(name_fo) ;

	  prompt("Enter name of file containing Wilson data for scaling:  " ) ; 

	  sscanf(terminp, "%s", name_fc) ;
 	  fprintf(fp_log, "Wilson data file name for scaling is %s\n", name_fc) ;

	  factor = scale_fo(name_fo, name_fc) ;

	  sprintf(message, "\n scaling fobs") ;
          printTwice(message) ;
	  scalefobs( Folist, exp( factor ) );

          extend_filename(pwd_name, out_filename, "_apo");

          sprintf(message, "\nWriting scaled and apodized file, %s", out_filename) ;
          printTwice(message) ;
          if (HDFOUT)
             write_HDFfobs(out_filename) ;
          else
             writeFolist(out_filename) ;
     }
     else{
	sprintf(name_fo, "%s_wil", pwd_name) ;
	make_selected_wilson(name_fo) ;

	prompt("Enter name of file containing Wilson data for scaling:  " ) ; 

	sscanf(terminp, "%s", name_fc) ;
	fprintf(fp_log, "Wilson data file name for scaling is %s\n", name_fc) ;

	factor = scale_fo(name_fo, name_fc) ;
        sprintf(message, "NOTE: use scaling factor %g in your Solve input file\n",
                factor) ;
     }
     return ;
}

  /*******************************************************
  Prepare a list of non-zero squared amplitudes (I),
  (sigma*I) squared values, and nusq values.
  We also need sigsq for apodization of sigmas.
  *******************************************************/

void	prepare_fo()
{
     real	fo_nusq(int) ;	/* in crysutil.c; uses Hlist etc	*/
     real	amp ;
     real	sig ;
     int	n ;

     /* Handle (h,k,l) triplets except for (000) */

     for (n = 1; n < Nfob; n++) {

          amp = *(Fap+n) ;
          sig = *(sigap+n) ;

          *(Int+n) = amp * amp ;
          *(ksq+n) = fo_nusq(n) ;
          *(sigsq+n) = sig * sig ;
          *(sigsqI+n) = 4. * amp * amp * sig * sig ;

     }			
}

  /*******************************************************
  Apply exponential factor to (h,k,l) triplets except for 
  (000). using the standard  PISQ * delta * eta * dr * dr,
  Round off very low amplitudes. Sigmas are apodized using 
  a, b and sratio array.  See analyze_sigsq().
  *******************************************************/

void	apply_apod_fo() 	
{	  
     real	fo_nusq(int) ;
     real	efac ;
     int	n ;
     real	Fmin ;
     real	thisF ;

     Fmin = amp0 * AMP_EPSI ;

     for (n = 1; n < Nfob; n++) {
	     
          efac = exp(-(delta * delfac * fo_nusq(n))) ;

          *(Fap+n) *= efac ;

          if (*(Fap+n) < Fmin) 
              *(Fap+n) = 0 ;

          if (useSig) {
	      thisF = *(Fap+n) ;
              *(sigap+n) = sqrt(*(sratio+n) * (asig + bsig*thisF*thisF)) ;
          } 
     }
}

  /*******************************************************
  Read the input fobs file using standard read function.
  Set up arrays in (h,k,l) reciprocal lattice triplets,
  *******************************************************/

int	apodfo_startup() 
{
     int	HDFdatalength(char *, char *) ;
     void	a_readInput() ;		/* read input parameters */
     void	fsetDsqLimit(int, char) ;
     int	readfobs(char *) ;	
     int        read_HDFfobs(char *) ;
     int	checkSigmas(real *, int) ;	

  /************************************************************
  Nlen is the file length of the input fobs data.
  ************************************************************/
     int	Nlen ;		
     int	n ;

     a_readInput() ;

     if (HDFIN)
          Nlen = HDFdatalength(sf_filename, "/fobs/Nfob") ;
     else
          Nlen = filelength(sf_filename) ;

     strcpy(message, "apodfo_startup") ;

     Fap = Folist =  (real *) e_malloc(Nlen*sizeof(real), message) ;
     sigap = Siglist = (real *) e_malloc(Nlen*sizeof(real), message) ;
     Hlist =  (float *) e_malloc(Nlen*sizeof(float), message) ;
     Klist =  (float *) e_malloc(Nlen*sizeof(float), message) ;
     Llist =  (float *) e_malloc(Nlen*sizeof(float), message) ;
     
     if (useSig)
        sratio   = (real *) e_malloc(Nlen*sizeof(real), message) ; 

  /************************************************************
  Note that Nfob is the actual number of fobs used (i.e., it
  excludes values at a resolution too high for the current 
  problem).
  ************************************************************/

     if (HDFIN)
          Nfob = read_HDFfobs(sf_filename) ;
     else
          Nfob = readfobs(sf_filename) ;
     
     useSig = checkSigmas(Siglist, Nfob) ;

     Int = (real *) e_malloc(Nfob*sizeof(real), message) ; 
     sigsqI = (real *) e_malloc(Nfob*sizeof(real), message) ; 
     ksq   = (real *) e_malloc(Nfob*sizeof(real), message) ; 
     sigsq = (real *) e_malloc(Nfob*sizeof(real), message) ; 

     if (useSig) 
       fprintf(fp_log, "Sigmas will be used for weighting\n") ;
     else
       fprintf(fp_log, "Sigmas will not be used for weighting\n") ;

     fsetDsqLimit(Nfob, TRUE) ;

     n = (dsq_limit + 0.5*binwidth)/ binwidth + 1 ;

     sprintf(message, "# of bins set to %d, average bin size is %g.", 
		n, binwidth) ;
     printTwice(message) ;

     return (n) ;
}

  /*******************************************************
  Determine coefficients (a, b) for fitting sigsq as a
  linear function of Int -- i.e. minimize
	Sum [sigsq - (a+b*Int)]^2.
  Letting f to be the above expression, we set partial 
  derivatives of f wrt a and b to 0, then solve.

  We set up sratio to be an array of sigsq / (a + b*Int)
  *******************************************************/

void	analyze_sigsq() 
{
    int	   n ;
    real thisI, thissq ;
    real SInt = 0, SIntsq = 0, SsigsqInt = 0, Ssigsq = 0 ;
    real denom, numa, numb ;

    for (n = 0; n < Nfob; n++) {
	 thisI = *(Int+n) ;
	 thissq = *(sigsq+n) ;
	 SInt += thisI ;
	 SIntsq += thisI * thisI ;
	 SsigsqInt += thisI * thissq ;
	 Ssigsq += thissq ;
    }
    denom = SInt * SInt - Nfob * SIntsq ;
    numa = SInt * SsigsqInt - SIntsq * Ssigsq ;
    numb = SInt * Ssigsq - SsigsqInt * Nfob ;

    asig = numa / denom ;
    bsig = numb / denom ;

    sprintf(message, 
             "Sigma apodization coefficients are %g and %g", asig, bsig);
    printTwice(message) ;

    if (asig < 0) {
	 printTwice("Your sigmas are very irregular;") ;
	 printTwice("Please check the fobs input or turn off the useSig flag!") ;
	 Fatal("Quitting...") ;
    }

    if (bsig < 0)  {
	printTwice(
    "Warning: your sigmas are not increasing with reflection intensities;") ;
        printTwice(message) ;
        sprintf(message, "Resetting bsig from %g to 0", bsig) ;
        printTwice(message) ;
        bsig = 0 ;
    }

    for (n = 1; n < Nfob; n++) 
	*(sratio+n) = *(sigsq+n) / (asig + bsig * *(Int+n)) ; 
}
